function OUT = TUSIM_DGP(SEED,n,k,s,r2,alter,shuffle,depend,flat)
% Generate parameters for DGP for simulation.
% Last edit: 3/28/17 CBH
% OUT = TUSIM_DGP(SEED,n,k,s,scalecoef,alter,shuffle,depend)
% SEED - seed for random number generator
% n - simulation sample size
% k - number of x variables (excluding the intercept and treatment dummy)
% s - number of relevant x variables (excluding the intercept and treatment
% dummy).  Actual number of relevant variables is 4*ceil(s/4)
% r2 - target R^2 for regression of y on x
% alter - pass 1 for all coefficients the same sign and -1 for alternating
% coefficients
% trunc - k x 1 vector of truncation points
% alter - 1 for same coefficients, -1 to have alternating signs
% shuffle - 1 to have shuffle elements of covariance matrix, 0 to have
% fixed
% depend - 1 gives linearly declining correlations, 0 gives geometrically
% declining correlations
% flat - one gives a flt signal, 0 gives a signal evenly spread between
% coefficients proportional to 1/sqrt(s) and 1/sqrt(n) 
% Simulation model:
% y = a0 + x'*a1 + b0*d + (x.*d)'*b1 + e

rng(SEED);
OUT.SEED = SEED;

OUT.trunc = 1.28*rand(k,1);  % truncation point for putting point masses in x 
if shuffle == 1
    OUT.shuffle = randperm(k);   % permutation for covariance matrix
else
    OUT.shuffle = (1:k)';
end

OUT.alter = alter;
OUT.flat = flat;

OUT.n = n;  % sample size
OUT.k = k;  % raw number of x variables
OUT.s = 4*ceil(s/4);  % raw number of relevant variables
OUT.p = 2*k;  % total number of x variables
if depend == 1
    OUT.OMEGA = toeplitz(1-(0:(1/(2*k)):(1/2)-(1/(2*k))));
else
    OUT.OMEGA = toeplitz(.8.^(0:k-1));
end
OUT.S = chol(OUT.OMEGA(OUT.shuffle,OUT.shuffle));  % covariance matrix between latent controls
OUT.prop = .5;  % Treatment probability

% Model coefficients
if flat == 1
    OUT.a0noscale = 1;   % Overall intercept
    OUT.a1noscale = [ones(ceil(s/2),1);zeros(k-ceil(s/2),1)].*(((alter).^(0:(k-1)))');  %Coefficients on controls in interaction term
    OUT.b0noscale = 1;  % Coefficient on treatment dummy
    OUT.b1noscale = [ones(ceil(s/2),1);zeros(k-ceil(s/2),1)].*(((alter).^(0:(k-1)))');  % Coefficients on controls in leading term
    OUT.coefnoscale = [OUT.a0noscale;OUT.a1noscale;OUT.b0noscale;OUT.b1noscale];
else
    OUT.a0noscale = (1/sqrt(s));   % Overall intercept
    OUT.a1noscale = ((2/sqrt(s))*[ones(ceil(s/4),1);ones(ceil(s/4),1)/sqrt(n);zeros(k-2*ceil(s/4),1)]).*(((alter).^(0:(k-1)))');  %Coefficients on controls in interaction term
    OUT.b0noscale = (.5/sqrt(s));  % Coefficient on treatment dummy
    OUT.b1noscale = (4/sqrt(s)*[ones(ceil(s/4),1)/sqrt(n);ones(ceil(s/4),1);zeros(k-2*ceil(s/4),1)]).*(((alter).^(0:(k-1)))');  % Coefficients on controls in leading term
    OUT.coefnoscale = [OUT.a0noscale;OUT.a1noscale;OUT.b0noscale;OUT.b1noscale];
end

% Find scale of coefficients to produce desired R^2
nRep = 10000;
varz = 0;
for jj = 1:nRep
    x = randn(OUT.n,OUT.k)*OUT.S;
    X = (x - ones(OUT.n,1)*OUT.trunc').*(x > ones(OUT.n,1)*OUT.trunc');
    d = rand(OUT.n,1) > OUT.prop;
    w = X.*(d*ones(1,OUT.k));
    z = [ones(OUT.n,1),X,d,w];
    varz = varz + cov(z)/nRep;
end

R2diff = @(c) ((c*OUT.coefnoscale)'*varz*(c*OUT.coefnoscale))/(1 + (c*OUT.coefnoscale)'*varz*(c*OUT.coefnoscale))-r2;
scalecoef = fsolve(R2diff,1);

% Define scaled coefficients
OUT.a0 = scalecoef*OUT.a0noscale;  % Overall intercept
OUT.a1 = scalecoef*OUT.a1noscale;  %Coefficients on controls in interaction term
OUT.b0 = scalecoef*OUT.b0noscale;  % Coefficient on treatment dummy
OUT.b1 = scalecoef*OUT.b1noscale;  % Coefficients on controls in leading term
OUT.coef = [OUT.a0;OUT.a1;OUT.b0;OUT.b1];
OUT.TRUE = OUT.coef ~= 0;
OUT.scalecoef = scalecoef;

% Individual treatment effect
OUT.x0 = [0;zeros(k,1);1;.5*ones(k,1)];  % Point of interest for individual specific treatment effect
OUT.te0 = OUT.x0'*OUT.coef;  % True value for treatment effect for individual of interest

% Calculate expected profit differential of strategy where mail to everyone with b0+x'*b1 > c0
profit = zeros(10000000,1);
meany = zeros(10000000,1);
for ii = 1:10000000
    if ~mod(ii,1000000)
        disp(ii)
    end
    x = (randn(1,k)*OUT.S)';
    X = (x - OUT.trunc).*(x > OUT.trunc);
    d = rand > OUT.prop;
    w = d*X;
    meany(ii) = OUT.a0 + X'*OUT.a1 + d*OUT.b0 + w'*OUT.b1;
    profit(ii) = (OUT.b0+X'*OUT.b1)*(OUT.b0+X'*OUT.b1 > 0);
end
OUT.pi0 = mean(profit);
OUT.scaleNorm = mean((meany.^2));


